#ifndef THERMAL_FI_OPERATOR_H
#define THERMAL_FI_OPERATOR_H

#include BS_FORCE_PLUGIN_IMPORT ()
#include "fi_operator.h"
#include BS_STOP_PLUGIN_IMPORT ()
#include "thermal_calc_model.h"
#include "thermal_data.h"
#include "boundary_building.h"

namespace blue_sky{

	template <typename strategy_t, bool is_w, bool is_g, bool is_o>
	struct thermal_fi_operator_impl : fi_operator_impl<strategy_t, is_w, is_g, is_o>
	{

		typedef fi_operator_impl<strategy_t, is_w, is_g, is_o>  base_t;
		typedef typename strategy_t::index_t                    index_t;
		typedef typename strategy_t::rhs_item_t                 rhs_item_t;
		typedef typename strategy_t::item_t                     item_t;
		typedef thermal_calc_model<strategy_t>				    thermal_calc_model_t;
		typedef smart_ptr<thermal_calc_model_t,true>		    sp_thermal_calc_model_t;

		typedef calc_model <strategy_t>                         calc_model_t;
		typedef smart_ptr<calc_model_t,true>					sp_calc_model_t;
		typedef typename calc_model_t::data_t                   data_t;
		typedef typename thermal_calc_model_t::data_t           thermal_data_t;

		typedef typename thermal_calc_model_t::data_array_t     data_array_t;
		typedef smart_ptr<build_boundary<strategy_t>, true>		sp_boundary_builder;

		typedef typename strategy_t::csr_matrix_t           csr_matrix_t;
		typedef typename csr_matrix_t::i_vector_type		i_vector_type;

		thermal_fi_operator_impl (sp_thermal_calc_model_t &thermal_calc_model, sp_calc_model_t &calc_model, sp_reservoir_t &reservoir, 
			const sp_mesh_iface_t &mesh, sp_jacobian_t &jacobian, sp_jmatrix_t &jmatrix, data_dimens dimens, sp_boundary_builder boundary_builder):
		fi_operator_impl(calc_model, reservoir, mesh, jacobian, jmatrix),
		thermal_data_(thermal_calc_model->thermal_data), temperature_(thermal_calc_model->temperature),
		ss_diag_(jmatrix_->get_ss_diagonal ()),
		dimens_(dimens),
		boundary_builder_(boundary_builder)
		{			
		}

		void fi_operator_cell (const item_t dt, index_t i, item_t &total_volume){
			base_t::fi_operator_cell(dt, i, total_volume);

			thermal_data_t &data_i = thermal_data_[i];
			data_i.temperature = temperature_[i];
		}

		/**
		* \brief  Fills Ass part for cell i
		* \param  i Index of cell
		* \param  data_i
		* \param  ss_block
		* */
		void fill_ss_block (index_t i, const data_t &data_i, rhs_item_t *ss_block)
		{
			base_t::fill_ss_block(i, data_i, ss_block);
			index_t shift = n_sec_vars - 1;
			ss_block[shift] = volume_[i];
		}

		void
		fill_sec_rhs (index_t i, const data_t &data_i, rhs_item_t *s_rhs_block)
		{
			base_t::fill_sec_rhs(i, data_i, s_rhs_block);
			index_t shift = n_sec_vars - 1;
			sp_thermal_calc_model_t thermal_calc_model_(calc_model_, bs_dynamic_cast());
			s_rhs_block[shift] = -volume_[i] * (thermal_data_[i].temperature - thermal_calc_model_->thermal_old_data_.temperature[i]);
		}

		void save_prev_niter_vars (){
			base_t::save_prev_niter_vars();
			int is_line_search = calc_model_->ts_params->get_int (fi_params::SELECT_SOL_STEPS);
			if (is_line_search != 0)
			{
				sp_thermal_calc_model_t thermal_calc_model_(calc_model_, bs_dynamic_cast());
				thermal_calc_model_->thermal_prev_niter_data_.save (thermal_calc_model_);
			}
		}

		/**
		* \brief  Eliminates cell
		* \param  local_data_
		* \param  jac_block
		* \param  rhs_block
		* \param  sp_block
		* \param  s_rhs_block
		* */
		void eliminate_cell (const local_data &local_data_, rhs_item_t *jac_block, rhs_item_t *rhs_block, rhs_item_t *sp_block, rhs_item_t *s_rhs_block, rhs_item_t *ss_block)
		{
			// Schur complemet 1) App-ApsDssDsp and 2) Bp-ApsDssBs (Dss=Ass^(-1); Dsp=Asp)
			// 1)
			m_minus_mmm_prod(n_sec_vars, n_phases, jac_block, local_data_.ps_block, sp_block, ss_block);
			// 2)
			m_minus_vv_prod2(n_sec_vars,n_phases,rhs_block, local_data_.ps_block, s_rhs_block, ss_block);
		}

		inline void m_minus_mmm_prod(index_t n_sec_vars, index_t n_phases, rhs_item_t *jac_block, shared_vector <item_t> ps_block, rhs_item_t *sp_block, rhs_item_t *ss_block){
			for (int i = 0; i < n_phases * n_phases; i++)                            
			{                                
				int sum = 0;                                              
				int shift_C = i % n_phases;                                   
				int shift_B = i / (n_phases) * (n_sec_vars);                     
				for(int j = 0; j < (n_sec_vars); j++)                    
				{                         
					sum += ps_block[shift_B + j] * sp_block[shift_C + j * n_phases] / ss_block[j];       
				}                                                       
				jac_block[i] -= sum;                                           
			} 
		}

		inline void m_minus_vv_prod2(index_t n_sec_vars, index_t n_phases, rhs_item_t *rhs_block, shared_vector <item_t> ps_block, rhs_item_t *s_rhs_block, rhs_item_t *ss_block){
			for(int i = 0; i < n_phases; i++)                             
			{                            
				int sum = 0;                       
				for(int j = 0; j < n_sec_vars; j++)                   
				{                  
					sum += ps_block[i  * n_sec_vars + j] * s_rhs_block[j] / ss_block[j];          
				}                                                     
				rhs_block[i] -= sum;                                        
			}
		}

		bool
		check_norm (index_t istart)
		{
			item_array_t &sec_sol = jmatrix_->get_sec_solution();
			index_t shift = n_sec_vars -1;
			item_t sec_sol_err = 0;
			
			typename i_vector_type::iterator it;
			i_vector_type non_boundary_cells = boundary_builder_->get_non_boundary_cells();
			for (it = non_boundary_cells.begin(); it < non_boundary_cells.end(); ++it){
				index_t cell = *it;
				sec_sol_err += fabs(sec_sol[cell * n_sec_vars + shift]);
			}

			bool result = sec_sol_err < 1E-4; 
			return base_t::check_norm(istart) && result;
		}

		/**
		* \brief  Calculates additional previous variables
		**/	
		void calc_prev_variables(){
			for (index_t i = 0; i < n_cells_; ++i)
			{
				thermal_data_[i].prev_temperature = thermal_data_[i].temperature;
			}
		}

		/**
		* \brief  Calculates flux part of Jacobian
		* \param  dt
		* \return True
		* */
		bool
		block_connections_mpfa (const item_t &dt);

		rhs_item_array_t              &ss_diag_;		

		data_array_t                  &thermal_data_;
		const item_array_t            &temperature_;	

		sp_boundary_builder			  boundary_builder_;

		data_dimens dimens_; 
	};
}

#include "thermal_fi_operator_block_connections_mpfa_2.h"

#endif //THERMAL_FI_OPERATOR_H